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STAR ADAPTATION FOR TWO ALGORITHMS 


USED ON SERIAL COMPUTERS 

By Lona M. Howser and Jules J. Lambiotte, Jr. 

Langley Research Center 

SUMMARY 

Two representative algorithms used on a serial computer and presently executed 
on the Control Data Corporation 6000 computer were adapted to execute efficiently on the 
Control Data STAR- 100 computer. Gaussian elimination for the solution of simultaneous 
linear equations and the Gauss-Legendre quadrature formula for the approximation of an 
integral are the two algorithms discussed. This paper describes how the programs were 
adapted for STAR and why these adaptations were necessary to obtain an efficient STAR 
program. Some points to consider when adapting an algorithm for STAR are discussed. 
Program listings of the 6000 version coded in 6000 FORTRAN, the adapted STAR version 
coded in 6000 FORTRAN, and the STAR version coded in STAR FORTRAN are presented 
in the appendices. 


INTRODUCTION 

Many algorithms which are presently used on the Control Data Corporation 6000 
computer and executed in a serial mode are suitable for the Control Data STAR-100 com- 
puter. However, if these algorithms were converted line-by-line to the STAR coding, it is 
not likely that they would make an efficient STAR program and might actually produce code 
which would run inefficiently on a vector computer such as STAR. The 6000 code will need 
to be adapted for STAR to produce an efficient STAR code. This paper discusses two algo- 
rithms of this nature. A comparison of the 6000 coding and STAR coding of the identical 
algorithm is made for two different algorithms: one for the solution of simultaneous 
linear equations using Gaussian elimination with partial pivoting and the other for numer- 
ical evaluation of an integral using the Gauss-Legendre quadrature formula. The paper 
discusses how the 6000 program was adapted for STAR, the reasons for the adaptations, 
and some 'factors which should be considered when adapting a program. FORTRAN 
codings of the algorithms are presented in the appendices. The STAR codings use the 
FORTRAN language defined in reference 1. 


/ 



AIDS FOR ADAPTING AN ALGORITHM 


When beginning to adapt an algorithm for STAR, it is not enough to look at just seg- 
ments in the 6000 coding, but it is necessary to look at the entire algorithm to get the 
total picture. Some questions to pose are: What is the final result?, What is needed at 
various steps in the algorithm ?, and What is computed independent of other steps and 
what is repeated? 

It will be helpful to review a few definitions and terms which are important to 
remember when formulating a program for the STAR computer. (For a comprehensive 
discussion of STAR architecture and hardware instructions, see ref. 2.) 

(1) A vector is a set of elements stored in contiguous locations in memory. 

(Array and vector are used interchangeably in this discussion.) 

(2) Vector timing is the time required by the central processing unit to process a 
vector. It is obtained by the equation T = S + Z/c, where 

T time in clocks (1 clock represents 40 nanoseconds) 

S startup time, different for each vector instruction or macro 

I length of vector (number of elements in the vector) 

c constant depending on the type of instruction and whether the 

vector elements are each stored in 32- or 64-bit words 

(3) A page is a block of storage which contains data or instructions. A program is 
made up of one or more pages. (See ref. 3.) 

(4) A program's working set is the smallest set of pages which must be in central 
memory for the program to operate efficiently, 

(5) A page fault occurs when a program references a page which is not contained in 
central memory; that is, it is not in the program's present working set. 

(6) Paging is the process of bringing a page into central memory or releasing a 

page. 

When adapting the 6000 code to the STAR code, the following factors are important: 

(1) Use vector instructions. This requirement may mean reordering steps in an 
algorithm so that elements in contiguous locations can be operated on or it may mean 
rearranging storage. 

(2) Use long vectors. If a choice is available whether to use many short vectors 
or a few long vectors, use the long vectors unless the overhead to create the long vectors 
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is too great. This procedure reduces the effect of the startup time associated with each 
vector instruction. 

(3) Avoid or use sparingly coding which will generate costly vector instructions 
and macros, that is, costly as compared with some of the faster vector instructions. 

Refer to the most current timings available. Examples of such instructions and macros 
are divide, transmit indexed list, and dot product. 

(4) Avoid unnecessary paging problems by creating a reasonable working set for 
the program. When working with a particular array, perform all the operations possible 
with this array before working with another array. This factor will be more critical 
with long arrays, but it is a good habit to form and should help reduce page faults. 

(5) Investigate the feasibility of creating more answers than are really needed. 
Because of the high result rate of vector instructions, it may be advantageous to use an 
approach which generates a larger number of results than are needed in order to avoid 
scalar computation. This idea, however, should be used cautiously. (See the discussion 
in ref. 4 on parallel algorithms for tri-diagonal equation solvers.) 

ADAPTATION OF AN ALGORITHM FOR THE SOLUTION OF 
SIMULTANEOUS LINEAR EQUATIONS 

A Langley Research Center 6000 library subroutine GELIM (see listing in appen- 
dix A) uses Gaussian elimination with partial pivoting to obtain the solution of the set of 
simultaneous linear equations, AX = B, where A is the square matrix of coefficients 
of order n, X is a vector of unknowns of length n, and B is a constant matrix of 
order n X r where r is the number of right-hand sides. The matrix A is factored 
into a lower unit triangular matrix and an upper triangular matrix. (For numerical 
details of the algorithm, see any numerical analysis text, such as ref. 5.) The subroutine 
also contains an option for the evaluation of the determinant of matrix A. 

The version of GELIM on the Langley Research Center 6000 library performs 
Gaussian elimination as generally defined by operating on one row of A at a time. At 
the kth step, column k is searched for the largest element; then row k and the row which 
contained the largest element are interchanged. The pivot element is then used to obtain 
zeroes in all positions in its column below the diagonal. This procedure requires multi- 
plying, a row of A by the appropriate scalar and subtracting this product from another 
row of A. 

Since these row modifications are the usual way the steps in Gaussian elimination 
are thought of being performed, it would seem normal that for STAR the matrix would be 
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stored by rows. This method of storage means that elements in one row of a matrix 
would be stored in contiguous locations so a row of the matrix could be a vector and the 
steps would be performed by using vector instructions. 

In the present version of STAR FORTRAN, two-dimensional arrays are stored 
columnwise and an optional storage arrangement by rows is not available. This column 
storage means that elements in one row of matrix A are not stored in contiguous loca- 
tions and modifying the matrix a row at a time would mean that few vector instructions 
could be used. This fact makes it desirable to see whether Gaussian elimination can be 
performed by modifying matrix A by using columns so that vector instructions can be 
used. As will be shown below, Gaussian elimination can be performed by operating on 
columns of the matrix. 

An option to store matrices by rows may be available in a later version of 
STAR FORTRAN, but when the 6000 code was modified to perform Gaussian elimination 
by columns, many advantages for the column storage over the row storage appeared; 
therefore, column storage is recommended. The following section will show how 
Gaussian elimination can be performed by using vector instructions when the matrix is 
stored by columns and will identify these advantages. The section entitled "Row Storage" 
shows the sequence of steps performed in the row storage which is identical to the pres- 
ent 6000 algorithm. 


Column Storage 

Gaussian elimination can be performed by modifying one column of the matrix at a 
time. This is done by a reordering of the operations from the usual row operations. 
Accomplishing the triangularization of matrix A by performing the work on columns 
makes efficient use of STAR and does the identical arithmetic normally done when per- 
forming the work on rows. 

The kth step of the triangularization can be performed as shown below, where n 
is the number of equations and r is the number of right-hand sides. All references to 
the kth column refer to column entries below the diagonal. 

(1) Divide the kth column of A by the aj^ element and store in the kth column. 
This is a vector divided by a scalar: 


_ ^ik 
ik atk 


(i = k + 1, . . ., n) 


(2) Multiply the kth column by the aj^j element and subtract this result from the 
jth column. This is a vector multiplied by a scalar and then a vector subtract: 
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~ " ^ik%j 


(i = k + 1, . . n) 


(3) Repeat the sequence of vector instructions in step (2) for the columns of 

A (j = k + 1, . . n). 

(4) Multiply the kth column of A by the element of B and subtract this 

result from the jth column of B. This is a vector multiplied by a scalar and then a 
vector subtract: 

bij = bij - aikbkj (i = k + 1, . . n) 

(5) Repeat the vector instructions in step (4) for the remaining right-hand sides 

of B (j = 1, 2, . . r). 

(6) Repeat steps (1) to (5) until the triangular ization of matrix A is complete 
(k = 1, 2, . , n - 1). 

Often a subroutine for the solution of simultaneous equations requires the user to 
append the right-hand sides to the original matrix. The column storage arrangement 
allows the right-hand sides to be done as separate vector instructions without having to 
append the right-hand sides to the original matrix. This arrangement is less cumber- 
some for the user since the right-hand side can be a separate array. 

GELIM uses partial pivoting which means that rows will need to be interchanged at 
times. At the kth step, column k is searched for the largest element and the row con- 
taining the largest element and the kth row are interchanged. The column search to find 
the largest element can easily make use of vector instructions, but the row interchange 
presents a problem. Elements of rows will need to be interchanged and none of them 
will be stored in contiguous locations. Not only are they not contiguously stored, but for 
a very large matrix a column could use one or more pages. For a large matrix it is 
unlikely that a program will be allowed a working set large enough to contain the entire 
matrix. This situation would mean that the row elements to be interchanged would be on 
separate pages and would have to be brought into and then out of core only to reference 
two elements on the page. Then when the column modifications are performed, the same 
pages will need to be brought back into core again. This procedure would be extremely 
inefficient. 

A form of indexing could be set up to achieve row interchange, but it would need 
one of the more time-consuming instructions (transmit indexed list) and the referencing 
across page boundaries would still be present. No vector instruction can help perform 
the interchange efficiently; thus, a scalar interchange will be just as efficient. 
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The possible paging problem can be alleviated by not interchanging an entire row 
at one time as it is performed on the present 6000 version. At the kth step when the 
largest element in column k is found, interchange the two elements only in column k. 

When working with the jth column, both column k and column j are needed. Both of the 
elements which need interchanging are in the jth column which will be in the program's 
working set at that time. Therefore, before beginning the operations on the jth column, 
interchange the two elements. In this way, the page or pages containing that column will 
need to be brought into core only once per step of the algorithm. 

The subroutine also computes the value of the determinant of matrix A which is 
equal to the product of the elements of the diagonal of the triangular matrix. In the 
6000 version, the product is computed after the triangularization is completed. For the 
STAR, this computation presents a similar situation as the rows interchange; the diagonal 
elements may be on separate pages and are not stored contiguously. Therefore, vector 
instructions cannot be utilized. The product will be scalar multiplication, but after the 
kth step has been completed, the partial product can be formed by using the diagonal ele- 
ment of column k while column k is still in the program's working set. This grouping of 
row interchange, triangularization of the matrix, and evaluation of the determinant should 
create an efficient working set for the program with a minimum of paging. 

The remaining task of the subroutine is to perform the back substitution. Back 
substitution generally uses the dot product of a row and the solution vector. This method 
is used on the 6000 version, but presents problems for STAR column storage since it 
uses the costly dot product macro and references elements of one array by rows and 
references elements of the other array by columns. This problem can be eliminated by 
reordering the steps needed to perform back substitution and all the work can be per- 
formed on columns. 


The steps to find the kth unknown by back substitution by using vector instructions 
but not using the dot product are as follows: 

(1) A scalar divide is always necessary, b, - = This step obtains the kth 

^ ^kk 

unknown for the jth right-hand side and stores the unknown in the right-hand side vector. 

(2) Multiply column k of A by the kth unknown obtained in step (1) and subtract 
this result from the jth column of B. This is a vector multiplied by a scalar and then a 
vector subtract: 


^ij = bij - aik^kj {i = 1, 2, . . ., k - 1) 

(3) Steps (1) and (2) are repeated for all the right-hand sides (columns of B) 
(j = 1, 2, . . ., r). 
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(4) Steps (1), (2), and (3) are repeated for all the unknowns (k = n, . . 2, 1). 

When k = 1, step (2) is omitted. 


Row Storage 

The way of performing the steps in Gaussian elimination as commonly seen in texts 
is by operating with rows. The elements in one row would be stored contiguously and 
the triangularization of matrix A would make efficient use of vector instructions. 

To make the triangularization most efficient, the right-hand sides must be appended 
to the matrix A for row storage. The length of the longest vector used would be 
n + r - 1 where n is the number of unknowns and r is the number of right-hand 
sides. If a separate array is used for the right-hand side, two identical vector instruc- 
tions would be needed for each operation, one of length n and one of length r. This 
method would be inefficient because the startup times would be multiplied by a factor 
of 2 and if r is 1, it would mean using vector instructions of length 1. Appending the 
right-hand side is a disadvantage in that it is awkward for the user. 

Step k of the triangularization when the matrix is stored by rows is as follows: 

(1) Perform a scalar divide, and store in a^^. 

(2) Multiply row k of A by a^ and subtract from the ith row. This is a vector 
multiplied by a scalar and then a vector subtract: 

Hi = Hi - ^ik^kj (j = k + 1, . . ., n + r) 

(3) Repeat steps (1) and (2) for all rows (i = k + I, . . ., n). 

(4) Repeat steps (1), (2), and (3) for all columns until the triangularization is com- 
plete (k = 1, 2, . . ., n - 1). 

The row interchanges necessary for partial pivoting are very easy when the matrix 
is stored by rows, but the column search will be scalar operations. In addition, there 
will be no way to avoid possible paging problems for a large matrix during the column 
search. The elements on each row of the column may be on separate pages, but the 
search has to be completed before any row operations can begin. 

The determinant evaluation could be performed the same way as in the column 
storage, after a step in the Gaussian elimination is completed. The pages for a large 
matrix would not have to be brought in only to get one element for the evaluation of the 
determinant, but the information would be used while the page was still in memory. 

When performing the back substitution, if the matrix is stored by rows, there is no 
way of using the scheme devised to eliminate the need for the dot product macro. The 
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back substitution references a row of the matrix and a column of the right-hand side. To 
be able to use vector instructions, a vector would need to be created which would contain 
the elements in the column. This method would be expensive and inefficient, and it is 
likely that the vector code would be no better than the scalar code. Either would be 
inefficient here, because of the referencing of a row and a column. 

Table I summarizes the comparison of the two storage arrangements at the various 
steps in the algorithm. The column storage is more advantageous for STAR than the row 
storage arrangement. 


TABLE I.- COMPARISON OF ROWWISE AND COLUMNWISE 
STORAGE OF MATRIX 


Step 

Rowwise 

Columnwise 

Column search 

Scalar: no way to avoid 
possible paging problems 

Vector 

Row interchange 

Vector: easy 

Scalar: can avoid possible 
paging problems 

Triangularization 

Vector: usual way it is done 

Vector: with steps reordered 

Back substitution 

Dot product macro: costly, 
referencing row and col- 
umn, possibly not 
vectorizable 

Vectors: columns only, no 
dot product, more efficient 

Determinant evaluation 

Scalar 

Scalar 

Treatment of right-hand 
sides 

Must be appended to matrix 
for efficient vector use 

Can be separate variable and 
still use vectors efficiently 


Flow Chart 

Subroutine GELIM was adapted for use on the STAR computer by using the same 
numerical method that was used on the 6000 computer. The matrix is stored by columns. 
By reordering computational steps, vector instructions can be used and a reasonable 
working set for the program established. Figure 1 shows a flow chart of the 6000 and 
STAR versions of the algorithm. 
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Figure 1,- Flow chart of 6000 and STAR version of segments of Gaussian 
elimination. An asterisk denotes use of vector instructions. 






















60C0 Back Suustilution STAR 



Figure 1.- Concluded. 


10 










Appendix A contains coding of the versions of the algorithm. The 6000 version and 
the STAR version coded in 6000 FORTRAN and in STAR FORTRAN are included. The 
calling sequence of the subroutine remains the same and no additional storage was 
necessary. 

Additional Comment 

An additional interesting fact was noted when the STAR version of subroutine 
GELIM was executed on the 6000 computer. The STAR version ran faster than the 
6000 version; as a result, there was about a 10-percent decrease in execution time. 

ADAPTATION OF A NUMERICAL INTEGRATION ALGORITHM 

Subroutine GLEGEN (see listing in appendix B) uses the Gauss- Legendre quadrature 
formula to evaluate simultaneously an array of integrals 

pb 

dx (k = 1, 2, . . N) 

where N is the number of functions. This subroutine is an example of the manner in 
which vectors rather than single variables can be used on STAR. As in the preceding 
example, the sequence of instructions is not the same as in the 6000 version. The subrou- 
tine subdivides the integration interval (a, b) into NQ panels and the Gauss- Legendre 
quadrature formula is applied to each panel by using a 3- or 10-point formula. The 
resulting integral is 

NQ NP 

•k = ^ J 2 N) 

i=l 1=1 

where 

fj^ function 

NQ number of quadratures or panels 

NP number of points per quadrature 

p^ = + AXj^ (where is the lower limit of integration for quadrature j) 

r- weights for the point formula used 
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abscissas for the point formula used 


b - a 
NQ 


(For a complete discussion of the Gauss- Legendre quadrature formula, see any numeri- 
cal analysis text, such as ref. 5.) 

The 6000 version of the subroutine obtains the first point in a quadrature, evaluates 
the function at that point, multiplies the fxmction by the appropriate weight, and adds the 
product to a sum. It continues in the same manner until all the points have been evalu- 
ated and used in that quadrature; then the process is repeated for the remaining quad- 
ratures imtil the integral has been evaluated over the desired range. The sum is then 
multiplied by the appropriate delta to obtain the value of the integral. 

Looking at this sequence of instructions shows that much of this work is actually 
computed independently. The STAR version takes advantage of this. All the information 
required to compute the points is known so all the points for all the quadratures can be 
computed to form a long array. This array can be passed to the subroutine to evaluate 
the functions and a savings can be obtained by using vector instructions to evaluate the 
functions. This routine will now return an array of function values. All the weights are 
known so an array can be formed which contains the weights for all the quadratures. 

The remaining task is to multiply the weights by the functions and sum the results. This 
computation can use the dot product macro; even though it is costly compared with most 
vector instructions, it is less expensive than a vector multiply followed by scalar adds. 
This formulation requires the use of the dot product only once. A scalar multiply of the 
sum by the appropriate delta completes the integral evaluation. 

The weights and abscissas for each quadrature are identical; therefore, the length 
of these arrays in the 6000 version is equal to 10 for the 10-point formula. A version of 
the subroutine using short vectors could have been formed. The STAR version could 
have used vector instructions and computed results using one quadrature at a time with 
the same array of weights and array length as in the 6000 version. Computing the results 
one quadrature at a time means that the startup time associated with each vector instruc- 
tion and the startup times used in the function evaluation also would be multiplied by the 
number of quadratures used. Depending upon the nature of the function, this amount of 
time could be significant. 


Vector transmits could be used to create a long vector of weights which is really a 
repetition of the short vector. The long vector is 

q repetitions 


(ri, rg, . 




m’ 


ri, rg, 


•' ^m) 
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'where m depends upon the quadrature formula and q is the number of quadratures. 
The cost of setting up this vector seems to be minimal when compared with the multipli- 
cation of the startup times which would have prevailed in the short-vector version. 

A problem of this type (excluding the function evaluation) would probably not pro- 
duce a significant paging problem. All the arrays are one dimensional and will most 
likely not be extremely long. However, consideration should still be given to the working 
set which is created. Currently, there is a tendency to initialize variables at the 
beginning of the program. The vector of weights could be created at the beginning of 
the program, but they are not used until after the functions are evaluated; therefore, the 
vector is not created until just before it is used. Forming the vector of weights and then 
immediately using them results in less paging; therefore, a better working set is created 
than would have been created by initializing the vector. 

Subroutine GLEGEN was adapted for the STAR computer by using the same numeri- 
cal method as was used on the 6000. An efficient STAR routine was obtained by using 
long vectors and reordering the sequence of instructions. Figure 2 shows a flow chart 
of the 6000 and STAR versions. 

Appendix B contains the coding of the two versions of the algorithm. The 6000 ver- 
sion and the STAR version in 6000 FORTRAN and STAR FORTRAN are included. Notice 
that the calling sequence of the STAR version reflects the use of long arrays in that more 
working space is needed. 


CONCLUDING REMARKS 

Two algorithms used on the Control Data Corporation 6000 computer were adapted 
for use on the STAR computer. This adaptation required a rethinking of the flow of the 
entire problem. Array variables are used where the 6000 code used a single-value 
variable and the steps in the algorithm are not computed in the same order for the two 
versions. This reordering of steps and changes in variable assignments allow vector 
instructions to be used and a reasonable working set can be established for the program; 
thereby efficient use of the STAR computer and some implied improvements for the 
Control Data Corporation 6000 computer are made. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va,, February 20, 1974. 
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Figure 2.- Flow chart of 6000 and STAR versions of Gauss -Legendre 
quadrature formula for numerical integration. An asterisk 
denotes use of vector instructions. 
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APPENDIX A 


FORTRAN CODING OF 6000 AND STAR VERSIONS OF ALGORITHM FOR THE 
SOLUTION OF A SET OF SIMULTANEOUS EQUATIONS 


6000 Version in 6000 FORTRAN 




?^UR WOUT I NE GEL I M ( A,N,B t iSiKHS f MAXhf IPI VOT , 1 OP » U t: T lRM ^ J SC A L t ) 

GELIMOU i 



0 TMe\S I ON A ( MAXN t f'^AXN ) • B ( MA XN • i>Jp HS > ♦ I PI VCT < M AXN ) 


G i: L 1 0 U 2 

- 


A(MAXN*MAXN) = SQUARE MATKiA JF C Ot F F I C I f NT S ( A IS DESTKOYFU) 

GELiN^03 



N = MU'IBPR OF «3WS AND COLUMNS IN A 


GEL IMGO^ 

- 


8(MAXN,NRHS) = •1ATRIX CF CONS TANT S< P £ PL AC ED 8Y SOLUTION 

MATkl X) 

GLLIM0O5 



N8US = NUMBER Cf- COLUMNS IN B 


GlL I hOOo 

- 


MAXN - maximum number OF ROWS AliD COLUMNS IN A 


GEL i Mu u 7 



IPIVOT(MAXN) = RFCURD OF =OW INTEkCHANOcS 


GlL IM uO b 

- 


IQP . T 0P = l . EVALUATE DETERMINANT, I0P = 0,SK.IP OETEkMINANT EVALUATION OjV 

- 


OETERM - GIVES VALUE OF uE TEk MI NA T ( Dli TE RM = ( 10 Y* 100 i ** I S C AL E*0c TL HM ) OiO 

- 


ISCALt = SCALt FACTCJR COMPUTtu BY SUriKOUT INh IN UtTbfRK 

LVALUATION 

GhLlMOil 



TO KERP ntTFRM WITHIN THE FLOATING POINT wGRU 

SIZE (jF the 0L2 



CO'^PUTEH 


GEL IMuii 



SIGN=l. 0 


GELIM020 



DO 35 I=lfM 


btL i MOi 0 

35 


tPf VOT( I )=I 


GLL IMu^O 


5 

1 SCALE = G 


GLL 1 MO jO 


6 

R 1= 10* 0»«100 


G cL I Muo 0 


7 

R2:= 1.0/R 1 


GEL1M070 

iO 

OtTFRi^= 1 .0 


GC-L 1 MOd C 



NML=N-l 


GtL IM090 



IF(NMi.EQ*0) GO Tr, lOOO 


GELIMI jO 

JNO LAKGfc5T PEMAIMi'JG Tt.P'1 IM I-TH CoLUKN ^■0t< PIVOT 


GcL IMilO 



GO 200 I = KNMi 


GEL IMI20 



R1G=0. 


GEL IM130 



DO 50 K= I ♦N 


GELIMJ.-4U 



TFRW=aBS( A(K,I ) ) 


GEL IMl^U 



IF(TEP'^'RIG> 50. 50.30 


GELIMIOO 

30 


RIG=T£RM 


GEL 1 Mi 7 w 



L = K 


GEL IMLdO 

50 


CQNTINIIF 


GEL l^li9G 



IPIVPT( T 


GEL lM2uO 



!F(RIG)BO. 60.30 


GEL IM210 

50 


D£TEPM-0*0 


GEL lM^<^Ci 



RL-TUPN 


GtL IM225 

ao 


I F( t-L> 90.120.90 


GLL IM2J0 

POWS Of A ANIj B* SCT IPI VD T ( J i - K ♦ 


obL iM2^0 

QO 


S!CN=-S1 GN 


GEL I M250 



DO 1 00 J = 1 . M 





TEMPl=A ( I, J) 


GEL I M2 7 0 



AlIiJ)=AlL.J) 


GEL IM2BG 

1 00 


A (L, JI = TEMPI 


GEL 1M.290 



DO ml .l = l,N‘^HS 


GLl IM iJu 



TEMP 2=8 1 I t 


uELIMjIO 



fl ( t . J ) = B ( L • J ) 


GEL IM32U 

101 


BIL. J)=TEmP2 


GELlMiJO 

120 


CONTINUE 


GELIMB^O 

PIVCT IN MUI.TIOLY A ANU 8 HY PIVOT. 


GEL IM35G 



IPl = H-L 


GEL lM3uO 



DO 31 I T=IP1,N 


GEL 1 M3 7 u 



A(ii,ii = A{ii.n/A(i.i) 


GELIM330 



X3=A( II , I 1 


GEL 1M3V0 



DO 32 X=IPl ,N 


GELlM^uO 



APPENDIX A - Continued 


32 A(I I .K» = 4(I ! ,Kl-X3=i>A( I,K) 
on 33 K=1,NRHS 

33 B ( I I ,K I -H ( I I .K )- X3*R ( I . K) 

31 COMTI^llig 

2Q0 CPNTINUE 

rF( Jf'P.eO.OIGC 221 

the OBTEBM INAMT 

on 122 r=i.N 

OIVOT1 = A(I , 1 ) 

1035 lF(AflS(l)ETEPM|-R l)10 3u,1010tl010 
1010 OETER'^=DETERM/P 1 
I':CALF= ISCALEH 

lF-lAB5(0tT6PM|-P 1)1060 tlU20. 1020 
1020 OETEP.M=r)ETEBM/R 1 
I 5CALE= 1 5CALfc* 1 

00 TO 1060 

1 030 ! HI A6S( DhTERM) -P2)l040i 10A0»lu60 
1040 OeTfcRH=DE‘^E»M*R I 

1 SC Al E= I SCALfc-l 

I FI ABS( DETERM) -R2) 10 50, 10 50, 1060 
1050 0£TEBM=0ETERM*R1 
I Sr ALE- I SCAL E- 1 

1060 IF! ABSI PIVOT 1 l-RUlO-yO, 1070, 1070 
1070 PtVOTI*PIVOT(/Rl 
I SCALE=I Sr.AL fcfl 

TF(ABS(PIVDTI)-Rn320,i080,10o0 
loao PivoTi=pivori/Ri 
I SCALE* ISCALE+1 
(Sn TO 320 

109 0 IF( ABS( PIVOT n-R 2 12000 , 2000, 320 

20CC PIVOTI=PIVOT1«R1 
I SCALE* ISCALE-1 

I HI AilSI P I VOT n-P2 12010, 2010,320 
2010 IvnTl =P I vrST t*R I 
I SCALE* 1 SCAL =-l 


GtLlM 410 
GEL 1 M 420 
GEL1M430 
GEL IM<,40 
GEL 1M450 
GEL I H460 
Gfc L i M h / j 

GLL 

GEL1K5O0 
GtL IM5LC 
GtL IMb^O 
GhLIH^iC 
GELIH^'+U 
GEL IM5!ju 
GEL IMfitiO 

gelim;>/o 

GELlMSdO 
GELIM^^O 
GtLIMoJu 
GEL IMoiG 
GELIM620 
GELIM63C 
GELIMO40 
GLLIM650 
GEL IMO60 
GELlMofu 
GEL IH6dU 
GEL IH6VU 
GEL IMZGu 
GEL1M710 
GELIM720 
GELIM730 
GELIMZ^O 
GEL1M75U 
GELIM/oO 
GELIM7/0 


320 

DETERM^OETEKM^P I VOTI 

GELIM7aO 

122 

C0NTINU6 

GEL1H790 

r 


GEtIMBOO 

♦•♦♦♦♦PERFORM SACK SUBSTITUTION 

GEUMBIO 

221 

CONTINUE 

GEUH820 


DO 57 IC“1«NRHS 

GELIM630 


BIN, ICI*B(N, ICI/AIN, N) 

6ELIN840 

57 

CONTINUE 

GEL1M85CI 


00 61 KK«1,NM1 

GELIN860 


I*N-KK 

GEL1M870 


ii=r*l 

GELIM880 


DO 61 J*1,NRHS 

GELIM900 


SUM*8(T. J) 

GEL1M910 


on 62 K*I1,N 

GEL1M920 

62 

SUMxSUM-All ,K)*BIK,J) 

GELIH930 

61 

Btl, JI*SUM/AII ,I 1 

GEL1M940 


OETERM=DETERM^SIGN 

GELIM950 


RETURN 

GELIM960 

1000 

IF(AI1,11.EO.OI GO TO 60 

GELIM970 


00 1500 J=1,NRMS 

GELIM975 

1500 

B(l, J»*BI1,J|/A( 1,11 

GEI.1M980 


OETERM*A(l.l) 

GEL1M990 


RETURN 

GEUM99 5 


END 

GELHIOOO 
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APPENDIX A - Continued 


STAR Version in 6000 FORTRAN 

G1L60010 
GIL60020 
6IL60030 
G1L60040 
GIL60050 
G1L60060 
GIL60070 
GIL60080 

C lOP - I0P=lf EVALUATE DETERMINANT. lOPsO.SKIP DETERMINANT EVALUATIOG1L60090 
C OETERM - GIVES VALUE OF OETERMINATIOETERM=( lO**100J**ISCALE*OETERMGlL60iOO 
C I SCALE = SCALE FACTOR COMPUTED BV SUBROUTINE IN DETERM EVALUATION GIL60U0 
; TO KEEP DETERM WITHIN THE FLOATING POINT WORD SIZE OF THEGIL60120 


- 

COMPUTER 




GIL60130 


10*0 




GIL60L40 


S[GN*1«0 




GIL60150 


00 3 I=l*N 




GIL60160 

3 

IPIVOT( I 




GIL60170 

5 

ISCALE*0 




GIL60180 

6 

Rl=l0,0^*100 




G1160I90 

7 

R2*1.0/Rl 




G1L60200 

10 

0ETERM=l-0 




GIt60210 






GIL60220 


IF (NMl.EO.O) GO TO 600 



GIL602J0 


NPl=N4-l 




G1L60240 

LARGEST REMAINING TERM IN I-TH 

COLUMN 

fOR PIVOT 

GLI60250 


DO 200 I>1*NM1 




G1L60260 


BIG<0. 




G1L60270 


DO SO K>I.N 




GIL602aD 


TERM>ABS(AtK« I \ t 




G1L60290 


IF1TERM-BIG)50.50.30 




GIL60300 

30 

nlG=TERM 




G1L60310 


L*K 




G1L60320 

50 

CONTINUE 




GIL60330 


IPIVOTI I l»L 




GIL603^0 


(F<BIG)80t60t80 




GI16Q350 

&0 

0£TERM=0.0 




GIL60360 


RETURN 




GIL60370 

80 

1F< I-U90.120.90 




G1L60360 

:**♦** interchange aEMENTS 

IN COLUMN I ONLY 



90 

SIGN=-5IGN 




G1L60400 


TEMPI = All.n 




GILbO^lO 


AU. n>A(L.I ) 




G1L60420 


AIL. n*TEMPl 




GR60430 

120 

CONTINUE 




G1L60440 

;*****ST0R£ pivot IN Ad.UI. MULTIPLY A AND B BY 

PIVOT* 

G1L60450 






GIL60460 

-♦♦♦♦♦obtain COLUMN OF PIVOTS 





DO 128 K=IPl.N 




GIL60470 

128 

A(K. n>A<K. I >ZAI 1. 1! 




G1L60480 

-♦♦♦♦♦INTERCHANGE ELEMENTS 

ONLY IN COLUMN 

J 




DO 132 J=1.NRHS 




GR60490 


TEMP1= BII.J) 




G1L60500 


B(I« J)=B(L. Jl 




G1L60510 


8IL. J)=TEMPl 




GIL60320 


DO 130 K*IPl,N 




G1L60530 

130 

61K. JI=B(K.J1 - AIK, 

I»*BI I.Jl 



GIL60540 

132 

CONTINUE 




GIL60550 

-♦♦♦♦♦INTERCHANGE ELEMENTS 

ONLY IN COLUMN 

J 




DOIAO J»IP1.N 




GR60560 


TEMPI = AII.JI 




GIL60570 


AII.J) = AIL. J> 




GIL60580 


AIL.JUTEMPl 




GR6G590 


□0135 K=IP1,N 




GIL60600 


SUBROUTINE GEL IMIA.N.8. NRHS.MAXN. IP! VOT t lOP.DETERM.ISCALE) 
DIMENSION AIMAXN.MAXNI .BlMAXNtNRHSI tlPI VOTIMAXNI 
: AIMAXN.MAXNI » square matrix of COEFFICIENTSIA IS DESTROYED! 

C N = NUMBER OF ROMS AND COLUMNS IN A 

: BtMAXN.NRHS) « MATRIX OF CONSTANTSIREPL ACED BY SOLUTION MATRIX! 

C NRHS = NUMBER OF COLUMNS IN B 

C MAXN = MAXIMUM NUMBER OF ROMS AND COLUMNS IN A 

: IPIVOT(MAXN) - RECORD OF ROM INTERCHANGES 
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APPENDIX A - Continued 


135 A(K«J)=A<K*J)-A(Kfl}^AU«il 
140 CONTINUE 

IF ( IDP. 60.01 GO TO 200 

^♦♦■►tSCAL6 THE DETERMINANT 

PIVOTI = A(I.n 

1005 IF 4 ABS(OET 6 RH)-RXIL 030 f 1010* 1010 
lOlO 0 ETERM= 0 €TERM/RI 
ISCALE=I$CALE «-1 

IF< ABS 4 DETERHJ-R II 1060 . 1020 f 1020 
1020 D£TERM=OETERM/Rl 
ISCALE^ISCALEfrl 
GO TO 1060 

L 030 IF( ABS(DETERM)-R 2 )i 040 . 1040# 1060 
1040 06 TERM» 0 ETERM^R 1 
ISCALE=ISCALE-l 

IF( A 8 S( DET 6 RMI-R 2 I 1050 , 1050*1060 
1050 0 ETERM=DETERM»R 1 
I SCALED I SC At E “1 

lObO IF 4 ABSI PIVOT n-Rl 1 1090 # 1070.1070 
1070 PIVOTlcPIVOTI/Rl 
ISCALE^ISCALEfrl 

I F ( A 8 S 4 PIVOT I l-R 11320 . 1080# 1080 
1080 PIV 0 TI==PIV 0 TI/R 1 
ISCALE*ISCALE ^1 
GO TO 320 


GIL 60610 

G 1 L 60620 

G 1 L 60630 

GIL 60640 

G 1 L 60650 

GIL 60660 

G 1 L 6067 Q 

G1L60680 

G 1 L 60690 

GIL 60700 

6 IL 60710 

GIL 60720 

GIL 60730 

GIL 60740 

G 1 L 60750 

GIL 60760 

GIL 60770 

G 1 L 60780 

GIL 60790 

GIL 60800 

G 1 L 60810 

G 1 L 60820 

GIL 60830 

GIL 60840 

GIL 60850 

G 1 L 60860 

G 1 L 60870 


1090 I F4ABS4 PIVOT n-R2l2000. 2000# 320 
2000 PIVOTI=PIVOTI^Rl 
ISCALE=ISCALE-1 

IF4 ABS4 PIVOT n-R 2 12010.2010# 320 
2010 PIVQT1«PIV0TI^RI 
I SCAtE-I SCALE-1 

-♦♦♦♦♦OBTAIN PARTIAL PRODUCT FOR DETERMINANT 
320 OETERM»OETERM*PI VOTI 

IF (I.NE.NMII GO TO 200 
IF (lO.EO.Ll GO TO 200 
PIV0TI^A4N.N) 

ID*1 

GO TO 1005 
200 CONTINUE 

-♦♦♦♦♦PERFORM BACK SUBSTITUTION 
340 CONTINUE 

00 450 KK=1#N 
K=NP1-KK 

^♦♦♦♦♦SCALAR DIVIDE 

DO 430 J*1»NRHS 
B(K« J)«6(K« J I/A4K.K) 

IF (K.EO.U GO TO 430 
KMl=K-l 

: ♦♦♦♦♦COLUMN K BY UNKNOWN AND SUBTRACT 
DO 420 I=liKMl 

420 B<I.J» s BCIfJt - ACI.KI^ B4K«JI 
430 CONTINUE 
450 CONTINUE 

DETERM<=OETERM*SIGN 

RETURN 

bOO IF ( A( I , D.EO.O) GO TO 60 
DO 800 J=l#NRHS 
800 BU# J1«B4I.JI/A( 1. 11 
OETERM^^AU.ll 
RETURN 
END 


GIL 60880 

GIL 60890 

GIL 60900 

GIL 60910 

G 1 L 60920 

G 1 L 60930 

GIL 60940 

GIL 6 U 950 

GIL 60960 

GIL 60970 

GIL 60980 

GIL 60990 

G 1 L 61000 

GIL 61010 

G 1 L 61020 

GIL 61030 

GIL 6 I 040 

GIL 6 L 050 

G 1 L 6106 G 

GIL 61070 

G 1 L 61080 

GIL 61090 

GIL 61100 

GIL 61110 

GIL 6 U 20 

G 1 L 6 U 30 

GIL 6 U 40 

GIL 61 I 50 

G 1 L 61 I 60 

GIL 61 L 70 

GlL 6 il 80 

GIL 6 U 90 

G 1 L 61200 

GIL 61210 
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APPENDIX A - Continued 


STAR Version in STAR FORTRAN 

GLiSOOlO 
GLIS002G 
GLIS0030 
GL1S0040 
GLIS0050 
GLISU06G 
GLIS0070 
GLisooao 

lOP - lOPaltEVALUATE DETERMINANT. IOP=O.SKlP DETERMINANT E V ALUAT I OGLI S0090 
DETERM - GIVES VALUE OF OETERMINATI 0£TERM=( 10**100) *• I S CALE^UETERMGLI SOLDO 
ISCALE = scale factor COMPUTED BY SUBROUTINE IN OETERM EVALUATION GLISOilO 
TO KEEP OETERM WITHIN THE FLOATING POINT WORD SUE OF THEGLIS0120 


: COMPUTER 


GL1S0L30 

10=0 


GLIS0L40 

SIGN=1.0 


GLI S0150 

DO 3 I=ltN 


GLI S0160 

3 IPIVOT(I) = I 


GLiSOUO 

5 ISCALE = 0 


GL1SOI0O 

6 Rl=10. 0**100 


GLIS0190 

7 R2=l.0/Rl 


GLISU200 

10 0ETERM=1.0 


GHS0210 

NM1=N-1 


GLIS0220 

IF (NMl.EO.O) GO TO 600 


GLIS0230 

NPl=NU 


GLIS0240 

;*****FINO LARGEST REMAINING TERM IN I-TH 

COLUMN FOR PIVOT 

GLIS0250 

DO 200 I^liNMl 


GLIS0260 



GL1S02/0 

: NEED VECTOR ABSOLUTE • MAXIMUM BECAUSE 

YOU ARE SEARCHING COLUMN 

GLISO20O 



GLISU290 

BIG^O« 


GL1S03U0 

00 50 K=I,N 


GLIS03I0 

TERM*ABS ( A(Kt n ) 


GLIS0320 

IF<TERM-BIG) 50» 50.30 


GL1S0330 

30 BIG=T£RM 


GLIS0340 

L*K 


GLIS0350 

50 CONTINUE 


GLIS0360 

IPI VOTt I )^L 


GLIS0370 

IF(BlG)d0.60.80 


GtISO30O 

60 DETERM^O-0 


GLI S0390 

RETURN 


GLIS04UU 

00 IFi I-LI90.120.90 


GLIS04L0 

90 SIGN=-SIGN 


GLIS0430 

|;*****INT£RCHANG£ ELEMENTS IN COLUMN I ONLY 


TEMPI* All, I) 


GL1S0440 

All . l) = Aa.I) 


GLIS0450 

AIL, n = TEMPl 


GLI $0460 

120 CONTINUE 


GLIS0470 

[:****#ST0RE pivot in ACUJI. multiply a and B by PIVOT* 

GLISO40O 

IP1*I*1 


GLI S0490 

^♦♦♦♦♦OBTAIN COLUMN OF PIVOTS 

AIIPLIN.!)* A( IP1:N. II /A! 1,1) 


GLIS0500 

;*****INTERCHANGE ELEMENTS ONLY IN COLUMN 

J 


00131 J=l,NRHS 


GLIS05L0 

TEMPI* B(I,JJ 


GLIS0520 

BII,U) = B(L,J ) 


GLIS0530 

BIL, J)=TEMP1 


GLIS0540 

131 B(IP1:N,J) = B(1P1:N,JI 'A(1P1:N.I) 

* B(IfJ) 

GLIS055Q 

r ***** interchange ELEMENTS ONLY IN COLUMN 

J 


00134 J=IP1,N 


GLIS0560 

tempi* AI!,J) 


GLIS0570 

All , J)-A(Lt J) 


GLIS0580 

AIL, J)=TEMP1 


GL1S0590 


SUBROUTINE GEL I M < A.N. 0 1 NRHS.MAXN, IPI VOT . lOP, OET ERM. I SCA LE » 
dimension A(MAXN.MAXN) «B4 MAXN*NRHS).IPI VOTIMAXN) 

AIMAXN.HAXNI SQUARE MATRIX OF COEFF IC I ENTS ( A IS DESTROYED) 

N > NUMBER OF ROWS AND COLUMNS IN A 

BIMAXN.NRHS) = MATRIX OF CONSTANTS! RE PL ACEO BY SOLUTION MATRIX) 
NRHS = NUMBER OF COLUMNS IN B 

MAXN » MAXIMUM NUMBER OF ROWS AND COLUMNS IN A 
IPIVOT(MAXN) = RECORD OF ROW INTERCHANGES 
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APPENDIX A - Concluded 


134 A(IPl:N.Jt» AdPUNtJ) -A(IPl;N.n* ACI»J» 
IF tlOP.EO.O) GO TO 200 

^♦♦♦**scale the determinant 

PlVOTI = AII,n 

1005 IFl ABS( 0ETERM)-RlJ1030t lOlOtlOlO 
1010 OETERH*oeTERH/Rl 
I SCALE* I SCALE^l 

IF( A6S( OETERMJ-R 1)1060. 1020. 1020 
1020 DETERM=0ETERM/R1 
ISCALE=ISCALE*1 
GO TO 1060 

1030 IF( A8S(0£TERH)-R2)1040. 1040. 1060 
1040 0ETERH=DETERM*R1 
ISCALE=!SCALE-1 

IF( ABS( OETERM)-R2)1050. 1050.1060 
1050 DETERM»DETERM*R1 
ISCALE=!SCALE-1 

1060 IFl ABSl P I VOTU-R 1)1090. 1070.1070 
1070 PIV0TI=PIV0TI/R1 
ISCALE=ISCALE^1 

TFI ABS< PIVOT I) -R 1)32 0. 1080.1080 
1080 PIV0TI=PIV0T1/Rl 
ISCALE = ISCALEU 
GO TO 320 

1090 IFI ABSI PIVOTI)-R2)2000.2000.320 
2000 PIV0TI = *>1VQTI*R1 
ISCAL£=ISCALE-1 

[FIABSI PIVQTI)-R2)2010.2010.320 
2010 PIV0T1=PIV0TI*RI 
I SCALE* I SCALE- I 

;****«>0BTA1N partial PRODUCT FOR DETERMINANT 
320 0ETERM=DETERM*PI VOTl 

IF ( I.NE.NMl ) GO TO 200 
IF (ID. EO.l) GO TO 200 
PIVOTI=A(N.N) 

10*1 

GO TO 1005 
200 CONTINUE 

*:perform back substitution 

340 CONTINUE 

00 450 KK=1.N 
K=NPl-KK 

C>*****SCALAR DIVIDE 

00 430 J*1.NRHS 
6(K* J)*e(K.J)/A(K.Kl 
IF (K.EO.l) GO TO 430 
KM1=K-1 

;**«.**C0LUMN K BT UNKNOWN AND SUBTRACT 

B(1:kmi.J)» BI1:kmi.J) -AI1:KH1.K)* BIK.J) 
430 CONTINUE 
450 CONTINUE 

DETERM*OETERM*S I GN 
RETURN 

600 IF (A(l.l).EO.O) GO TO 60 

B(1.1:NRHS)* B(1.1:NRHS)/AI1.1) 
D£TERM=Atl,l) 

RETURN 

END 


GL1S0600 
GLIS 0610 
GlTS0620 
GLIS0630 
6L IS 0640 
GL1S0650 
GLIS0660 
GLIS0670 
GLI S0680 
GLIS069G 
GL1S0700 
GLISOYIO 
GLIS0720 
GLIS0730 
GLIS0740 
GLIS0750 
GL1S0760 
GLIS0770 
GLIS07ao 
GLIS0790 
GLl SOUOO 
GLIS081U 
GLIS0820 
GLIS0630 
6LIS0840 
GLisoeso 

GLISOeoO 

GLIS087n 

GL1S0880 

GLIS0890 

GLIS0900 

GLIS0910 

GLIS0920 

6LIS0930 

GL1S0940 

GL1S0950 

GL1S0960 

GLIS0970 

GLIS0980 

GLIS0990 

GLISIOUO 

GLiSlOlO 

GLIS1020 

GLIS1030 

GLIS1040 

GL1S1O50 

GLIS1060 

GLIS1070 

GLIS1080 

6LIS1090 

GLISllOO 

GLISlllO 

GLIS1120 

GLIS1130 

GLISII40 

GLIS1150 

GL1S1160 

GLIS1170 


20 



APPENDIX B 


FORTRAN CODING OF 6000 AND STAR VERSIONS OF ALGORITHM 
FOR NUMERICAL INTEGRATION EVALUATION 

6000 Version in 6000 FORTRAN 


SUBROUTINE GLEGEN( SUM* I CODE » A»B« FX« FQFX • NFX # N0«NPQ J 

Z* * 

Z* PURPOSE: 

Z* 

Z^ 

Z* use: 

z* 

PARAMETERS: 

Z* SUM 


GLEGN 
GLEGN 
GLEGN 
GLEGN 
GLEGN 
GLEGN 

• GLEGN 

CALL GLEGEN < SUH« 1C0D£« At B* FXt FOFX f NF X»NOf NPO » « GLEGN 

« 

♦ 

THE ARRAY FOR THE YALUEISI OF THE INTEGRALISJ ♦ 


TD COMPUTE THE INTEGRALSt Fill OF X * OX FROM ♦ 
A TO B USING THE GAUS S-LEGENDRE QUADRATURE ♦ 
FORMULA^ * 


1 
2 

3 

4 

5 

6 
7 

a 

GLEGN 9 
GLEGN iO 
GLEGN II 


c ♦ 









4« 

GLEGN 

12 

:♦ 

ICODE 

AN 

INTEGER TEST CODE 

RETURNED BY 

THE 

ROUTINE 

♦ 

GLEGN 

13 

w * 


s 

0* 

NORMAL 

RETURN 





GLEGN 

14 

" * 


s 

1* 

NFX NOT 

PROPERLY 

SPECIFIED 



* 

GLEGN 

15 

:: * 


- 

2« 

NPO NOT 

PROPERLY 

SPECIFIED 

<NOT 

3 OR 101 

4c 

GLEGN 

16 

z* 


3 

3* 

NO NOT 

PROPERLY 

SPECIFIED 



4i 

GLEGN 

17 

z* 









4E 

GLEGN 

18 

3* . 


THIS 

PARAMETER SHOULD 

BE TESTED 

UPON 

RETURN 

4c 

GLEGN 

19 


z* 

A 

z* 

:♦ a 

z* 

Z* FX 

z^ 

z * 

:♦ FOFX 

z* 

Z* NFX 

:* NO 

z * 

Z* NPQ 

^ » 

Z ♦ 

Z* REQUIRED ROUTINES: 

Z* 

Z^ SOURCE/IMPLEMENTER: 

:* 

Z* DATE released: 

:♦ 

LATEST revision: 


BY THE CALLING PROGRAM 


LOWER LIMIT OF INTEGRATION 
UPPER LIMIT OF INTEGRATION 


GLEGN 20 
GLEGN 21 
GLEGN 22 
GLEGN 23 
GLEGN 24 
GLEGN 2b 

THE NAME OF A USER SUPPLIED SUBROUTINE WITH ♦ GLEGN 26 
ARGUMENTS X AND FOFX TO EVALUATE THE FUNCTIONS* GLEGN 27 
IT MUST BE DECLARED AS EXTERNAL* * GLEGN 28 

* GLEGN 29 


4t 

* 

♦ 

* 

* 

♦ 


THE ARRAY TO STORE THE VALUES OF THE FUNCTIONS* 

4c 

AN INTEGERf THE NO. OF FUNCTIONS TO INTEGRATE ♦ 

* 

AN INTEGER* THE NUMBER OF QUADRATURES * 

♦ 

AN INTEGER* THE NO. OF POINTS PER QUADRATURE ♦ 
IT MUST BE 3 OR 10 * 

♦ 

NONE * 

» 

« 

♦ 


COMPUTER SCIENCES CORP./ G. W. HAIGLER 
NOV. 14*1972 


NOV. 15*1972 


GLEGN 30 
GLEGN 31 
GLEGN 32 
GLEGN 33 
GLEGN 34 
GLEGN 35 
GLEGN 36 
GLEGN 37 
GLEGN 38 
GLEGN 39 
GLEGN 40 
GLEGN 41 
GLEGN 42 
GLEGN 43 
GLEGN 44 
GLEGN 45 
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APPENDIX B - Continued 


«**«««**«>«***********« *««***4'*«** **«*«*« *****«**«* GLEGN <»6 
Z GLEGN 47 

DIMENSION U1 (3 > «U2U0> t RU3I«R2( LOI «U(I3I tR<13i t SUM! 1). FOFXd I GLEGN 4S 


GLEGN 49 


equivalence (UK LKUI in , (U2m*U(4M*(«U(l> •(!( 1> I ,IR2( U.RI4>I GLEGN SO 

GLEGN 51 

DATA U 1 /, I 12 70 1 6653792 59.. 5t. 887298334620742/, UZ/. 01 304673574 14 14»GL£GN 52 
1.06 7468 3 1665 550 7,. 16029 5215650488.. 283302302935376,. 42556283050918GLEGN 53 
25, . 57443 71694908 16.. 716697697064624. . 839 704784149512, . 9325316 83344GLECN 54 

3493.. 986 953264258586/, R I/. 2 7777777777/778,, 444444444444444, .2777 77GLECN 55 
4777777778/, R2/. 033335672154344, .074725674575291,. 109543181257991, GLEGN 56 

5. 134633359654998.. 147762112357377.. 147762112357377. . 13463335965499GLEGN 57 

68.. 109543181257991, .074725674575291, .033335672154344/ GLEGN 58 


I CODE = 0 

IFlNFX.LE.OnCOOE = 1 
IFINO.LE.OUCOOE * 3 
IFINP0.E0.3)G0 TO 5 
IFiNPQ.NE.lOlICOOE = 2 
5 IF< ICOOE.GT.OIRETURN 
J = 3 

IF (NPO ,E0. 3> J=0 

DEIT= ( B-A)/FLOAT(NO» 
on 10 I = l,NFX 
10 SUM(I> > 0.0 
DO 80 K=l,NO 
XI » K - 1 
FF = A fr XI»OELT 
DO 80 L=1.NP0 
UU=U(J*L »*OELTfFF 
CALL FX(UU,FOFX) 

FACT»R( J*-U 
DO 80 J8=1,NFX 

80 SUM(JB» »FOFX( J8)*FACTfSUM(J6» 
00 100 I 1^1, NFX 
100 SUMUI1 = SUM(I!»*0£LT 
RETURN 
END 


GLEGN 59 
GLEGN 60 
GLEGN 61 
GLEGN 62 
GLEGN 63 
GLEGN 64 
GLEGN 65 
GLEGN 66 
GLEGN 67 
GLEGN 68 
GLEGN 69 
GLEGN 70 
GLEGN 71 
GLEGN 72 
GLEGN 73 
GLEGN 74 
GLEGN 75 
GLEGN 76 
GLEGN 77 
GLEGN 78 
GLEGN 79 
GLEGN 80 
GLEGN 81 
GLEGN 82 
GLEGN 63 


STAR Version in 6000 FORTRAN 


z* 

z* 

z * 

c* 

z* 

c* 

z* 
z* 
c* 
c * 
r.* 
z * 
z* 
z* 
z* 
z* 
z* 
z* 

c* 


SUBROUTINE GLEGE NS< SUM , 1 CODE. A. 8 , FX, FOFX , NFX.NO.NPQ.UU. FF . UO. MAXI 

* 

PURPOSE: TO COMPUTE THE INTEGRALS, F(II OF X • DX FROM ♦ 

A TO B USING THE GAUSS-LEGENORE QUADRATURE * 
FORMULA. ♦ 

* 

USE: CALL GLEGEN ( SUM, I CODE, A, B* FX.FOFX. NF X, NO.NPOl * 

* 


PARAMETERS: 


SUM 


I CODE 


THE ARRAY FOR THE VALUE(S) OF THE INTEGRALISi 


AN INTEGER TEST CODE 
» 0, NORMAL RETURN 
- I, NFX NOT PROPERLY SPECIFIED 
* 2, NPQ NOT PROPERLY SPECIFIED I NOT 
» 3, NO NOT PROPERLY SPECIFIED 


RETURNED BY THE ROUTINE 


3 OR 101 


THIS PARAMETER SHOULD BE TESTED UPON RETURN 
BY THE CALLING PROGRAM 


IGLENSOOl 

GLENS002 

GLENSU03 

GLENS004 

GLENS005 

GLENS006 

GLENS0G7 

GLENS008 

GLENS009 

GLENSOlO 

GLENSOll 

GLENS012 

GLENS013 

GLENS014 

6LENS015 

6LENS016 

GLENS017 

GLENS018 

GLENS019 

GLENS020 

GLENS021 
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APPENDIX B - Continued 


6 


FX 


FOFX 


NFX 


NO 


NPO 


UU 

FF 

UD 

MAXI 


REOUIPEO ROUTINES; 


LOWER LIMIT OF INTEGRATION 
UPPER LIMIT OF INTEGRATION 


♦ 
* 
* 

* 

A USER SUPPLIED SUBROUTINE WITH ♦ 
AND FOFX TO EVALUATE THE FUNCTIONS* 
iL. ♦ 

♦ 

A TWO OIMENSIONEO ARRAY TO STORE VALUES OF THE 
FUNCTIONSt FOFXIHAXIfNFXJ IN CALLING PROGRAM 


THE NAME OF 
ARGUMENTS X 
IT MUST ee DECLARED 


AN INTEGERf THE NO. OF FUNCTIONS TO INTEGRATE 
AN INTEGER, THE NUMBER OF OUAORATURES 


AN 

II 


INTEGER, THE 
MUST 8E 3 DR 


NO. 

LO 


OF POINTS PER QUADRATURE 


ARRAY TO STORE ALL POINTS AT WHICH FUNCTION IS 
EVALUATED ,ALSO STORES WEIGHTS 

ARRAY OF WORKING SPACEt CONTAINS LOWER LIMIT OF 
INTEGRATION FOR EACH QUADRATURE 
ARRAY OF WORKING SPACE. ABSCISSAS* OELT 
MAXIMUM ROM DIMENSIONS OF FOFX IN CALLING PROG. 


NONE 


SOURCe/IMPLEMENTER: COMPUTER SCIENCES CORP-/ G. W. HAIGLER 


GLENS022 
GLENS023 
GLENS024 
GLENS025 
GLENS026 
GLENS027 
GLENS028 
GLENS029 
GLENS030 
GLENS031 
GLENS032 
GLENS033 
GLENS034 
GLENS035 
GLENS036 
GLENSD37 
GLENS038 
GLENS039 
GLENS040 
GLENS041 
GLENS042 
GLENS043 
GLENS044 
GLENS04S 
GLE6S046” 
GLE6S047 
GLb6SQ48 
GLE6S049 
GL E 6S05 0 
GLE6S051 
GLE6S0&;; 
GLE6S053 
GLE6S054 

^t***********1^*****************l^*^^**^^*****^^*i^***^^*<^*1^*^l**^^^^^*t^*4l<^**^^ GLE6S05 3 

GLE6S056 
GLE6S037 
GLE6S038 
GLE6S039 
Gt.E6S060 

: GLE6S061 

DATA Ul /. 1 12 70 16653792 59.. 5.. 88729833462 C742/.U2/. 01 309 673 5 74 14 19.GLE6S 062 

1.067468316655507 .. 16029 521 5850483. .283302 3029353 76 .. 425 5628305091 8GLE6S06 3 

25.. 57443 7 1694908 16.. 716697697064624. . 839704784 1495 12 .. 9 325316 833 44GLE6S064 

3493.. 986953264258586/,ft 1/, 2777777777777 78,. 4444444444 44 444, ,277 777GLE6S06 5 
47 77 7777 78 /.R 2/, 03333 567 2154344,. 074725674575291,. 109543181257991, GLE6S066 

5. 134633359654998 .. 147762112357377.. 147762112357377. . 1346333 5965499GLE6S06 7 

68.. 1095 431 812579 91.. 07472 5674575291,. 0333356 72 154344/ GLE6S06 8 

: GLE6S069 

ICOOE = 0 GLE6S070 

NT = TOTAL. NUMBER OF POINTS AT WHICH FUNCTION WILL BE EVALUATED GLE6S071 

NT=NO*NPO GLE6S072 

IF( NFX.LE.OIICOOE - 1 GLE6S073 

IFINO.LE.O) ICOOE - 3 GLE6S074 

IF(NP0.E0.3IG0 TO 5 6LE6S075 

IFINPO.NE.lonCODE = 2 GLE6S076 

5 IF( ICODB.GT.OIRETUIN GLE6S077 

J=3 GLE6S078 

IF INPO .EO. 31 J=0 GLE6S079 

OELT= I B-AI/FLOATINOI GLE6S080 

;*****COMPUTE ARRAY TO BE ADDED TO FIRST POINT IN EACH QUADRATURE 

00 20 K=l,NPO GLE6S081 

J1 =J*K GLE6S082 

20 UD(K»= UUH*OELT GLE6S063 


DATE RELEASED: 


LATEST REVISION! 


NOV. 14,1972 


NOV. 15,1972 


* 

* 

* 

* 

* 

* 

* 


DIMENSION UI (31 ,U2( 101 ,R1( 3I.R21 10) ,U(1 3) ,R( 131 
1 FOFX(MAXI,l » ,UUU ) ,UD( 1 ),FF(1 ) 


SUM( 1), 


equivalence (U1 ( 1).U( 1)1 ,(U2( 1I.U(4) I, (R1(1),R( D) .(R2( 1),R(4) > 
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APPENDIX B - Continued 


C*****COMPUTE FIRST POINT IN EACH OUAORATURE 
00 30 K=1.N0 
XI =K-1 

30 FFiK) =A ^ Xl^OELT 

: CREATE VECTOR OF ALL POINTS FOR ALL QUADRATURES 

Du 4-0 K— i»NO 
K1 =^(K“l i^NPO 
00 35 I = 1*NP0 
35 UUlKlfl » = FFI KH- UO ( I ) 

AO CONTINUE 

: EVALUATE FUNCTION AT ALL POINTS 

C NOTE CHANGE IN CALLING SEQUENCE FOR FUNCTION EVALUATION SUBROUTINE 
CALL FXSCUUtFOFXfMAXNNT > 

: CREATE VECTOR OF WEIGHTS 

on 50 K=UNQ 
K1 =IK-l »*NPQ 
DO 45 1= 1,NP0 
45 UU(KH-I ) > RIJi-I I 
50 CONTINUE 

OQ 80 I=l#NFX 
SUM( n*o,o 
on 70 K=l*NT 

70 SUM(n*= SUMUl *UU(K|4 FQFXIK,U 
SUH{ II = DELT ♦ SUMC II 
80 CONTINUE 
RETURN 
END 


GLE6S084 

GLE65085 

GLE6S086 

GLE6S087 

GLE6S08B 

GLE6SO09 

GLE6S090 

GLE6S091 

GLE6S092 

GLE6S093 

GLE6S094 

GLE6S095 

GLE6S096 

GLE6S097 

GLE6S098 

GLE6S099 

GLE6S100 

GLE6S10I 

GLE6S102 

GLE6S103 

GLE6S104 

GLE65I05 

GLE6SI06 

GLE6S107 

GLE6S108 

GLE6S109 


STAR Version in STAR FORTRAN 


Z** 

z * 

z * 

z * 

z* 
z * 
z ♦ 
z * 
z * 
c * 
z * 
z • 

z ♦ 

z^ 

z* 
z * 
z* 
z* 
z* 
z* 
z * 
z * 
c ♦ 
c * 

♦ 

c * 


SUBROUTINE GLEGENSiSUHt ICODE* A,B, FX, EOF Xf NFXtNOt NPO»UU^ FFt UDt MAXI I 

« 

PURPOSE: TO COMPUTE THE INTEGRALSt Fill OF X ♦ DX FROM * 

A TO 6 USING THE G AUSS-LEGENORE QUADRATURE * 
FORMULA, * 

* 

USE: CALL GLEGEN ( SUM, I CODE f At FX^f QFX , NF Xt NO* NPU I ♦ 

♦ 


PARAMETERS: 


SUM 


ICOOE 


FX 


FQFX 


THE ARRAY FOR THE VALUEiSI OF THE INTEGRALISI 

AN INTEGER TEST CODE RETURNED BY THE ROUTINE 
* Of NORMAL RETURN 
= It NFX NOT PROPERLY SPECIFIED 
- 2, NPO NOT PROPERLY SPECIFIED (NOT 3 OR 10) 
= 3t NO NOT PROPERLY SPECIFIED 


THIS PARAMETER 
BY THE CALLING 


SHOULD BE TESTED UPON RETURN 
PROGRAM 


LOWER LIMIT OF INTEGRATION 


UPPER LIMIT OF INTEGRATION 


THE NAME OF 
ARGUMENTS X 
IT MUST BE DECLARED AS EXTERNAL, 


A USER SUPPLIED SUBROUTINE WITH 
AND FOFX TO EVALUATE THE FUNCTIONS^ 

A TWO DIMENSIONED ARRAY TO STORE VALUES OF THE 
FUNCTIONS, FOFXiMAXItNFX) IN CALLING PROGRAM 


GLE6SOOI 

GLE6S002 

GLE6S003 

GLE6S004 

GLE6SU05 

GLE6S006 

GLE6S0D7 

GLE6S006 

GLE6S009 

GLE6S0I0 

GLE6S011 

GLE6S012 

GLE6S013 

GLE6S014 

GLE6S015 

GLE6S016 

GLE6S017 

GL£6S0i8 

GLE6S0L9 

GLE6S020 

GLE6S021 

GLE6S022 

GLE6S023 

GLE6S024 

GLE6S025 

GLE6S026 

GLE6S027 

GLE6S028 

GIE6S029 

GLE6S030 

GLE6S031 
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APPENDIX B — Continued 


c* 

z* 

z* 

z * 
:♦ 
z* 
z* 
z* 
z* 
t* 


NFX 

NO 

NPQ 

UU 

FF 

UD 

_ MAXI 

REOUIReO ROUTINES: 


♦ GLE6S032 

AN INTEGER# THE NO* OF FUNCTIONS TO INTEGRATE ^ GLE6S033 

♦ GLE6S034 

AN INTEGER# THE NUMBER OF QUADRATURES * GLE6S035 

♦ GLE6S036 

AN INTEGER# THE NO* OF POINTS PER QUADRATURE « GLE6S037 
IT MUST BE 3 DR 10 * GLE6S038 

ARRAY TO STORE ALL POINTS AT WHICH FUNCTION IS GLE6S039 
EVALUATED #ALSO STORES WEIGHTS ^ GLE6S040 

ARRAY OF WORKING SPACE# CONTAINS LOWER LIMIT OF GLE6S041 
INTEGRATION FOR EACH QUADRATURE , GLE6S042 

ARRAY OF WORKING SPACE# ABSCISSAS* OELT GLE6S043 

GLE6S044 

MAXIMUM ROW DIMENSIONS OF FOFX IN CALLING PROG. GLE6S04^ 

GLENS046 

none * GLENS047 


NOV. 14#1972 


NOV. 15#1972 


♦ GLENS048 

♦ GLENS049 

♦ ' GUENSOSO 

♦ GLENS051 

♦ GLENS0&2 

♦ GLENS053 
GLENS054 
GLENS055 
GLENS0^6 
GLENS057 
GLENS058 
GLENS0!>9 
GLENS060 

DATA U1 /.I 1270 16653792 59#. 5# #887 29833462 0742/ #U2/. 01304673574141 4# GLENS061 
1.06 7468316655507# . 160295215850488. .233302302935 376,# 42556283050918GLENS062 
25#. 574437169490016#. 716697697064624# . 839704784149512# .9 32531683344GLENS063 
3493# -98695326425 8586/# R I/. 2 77 77 77 777 7 77 78 #.444444444444 444# • 2777 77GLENS064 
47 777777 78 /#R2/. 033335672 154344#. 07472 5674575291 ##109543 I 812 5799 I# GLENS065 
5. 1346333 59654996 #.147762 11235 7377## 147762 112357377## I 3463335963499GLEN S066 


I* 

;* SDURCE/IMPLEMENTERS COMPUTER SCIENCES CORP./ G. W. HAIGLER 
:* 

:• DATE released: 

;* 

LATEST REVISION: 

; #««««««««** ************************** ********** 

DIMENSION UIB) .U2( LO) ,Rl(3>.R2U0I.U(13>*R(li) .SUM( U« 

L FOFXIMAXI.l >,UU(l »#UO(l I.FFIl ) 

EQUIVALENCE (U1 1 U «U ( I U t ( U2I 1 1 .UU) I . ( R1 ( 1 ) >R( 1 ) ) f ( R2( 1 ) • R ( 4) > 


68,. 1095 43181 25 7991.. 07472 5674575291 ..033335672154344/ 


ICODE » 0 

; NT = TOTAL NUMBER 
NT=NO*NPO 

IF(NFX.LE.O) ICODE = 1 
IFINO.LE.OnCOOE = 3 
IFtNP0.EO.3)G0 TO 5 
IFINPQ.NE.IOICODE = 2 
5 IFnCODE.GT.ORETURN 
J*3 

IF (NPQ .EO. 3) J=0 

DELT* i e-A) /FLOATING) 

Jl = J*l 
NJ= J «-.NPQ 

;*****COMPUTE ARRAY TO BE 
UD(l:NPG) = U(Jl:NJ) 
:*****COMPUTE FIRST POINT 
00 30 K=l,NO 
XI *K-l 

30 FF(K) =A ♦ XUDELT 
;*«4r«*cOMPUTE ARRAY OF ALL 
00 40 K^L.NO 
Kl =<K-1I*NP0 


OF POINTS AT WHICH FUNCTION WILL BE EVALUATED 


ADDED TO FIRST POINT IN EACH QUADRATURE 
DELT 

IN EACH QUADRATURE 


POINTS FOR ALL QUADRATURES 


GLENS067 

GLENS068 

GLENS069 

GLENS070 

GLENS071 

GLENSQ72 

GLENS073 

GLENS074 

GLENS075 

GLENS076 

GLENS077 

GLENS078 

GLENS079 

6LENS060 

GLENS081 

GLENS082 

GLENS083 

GLENS084 

GLENS085 

GLENS086 

GLENS087 

GLENS083 

GLENS089 

GLENS090 


NK- K ♦ NPQ 

UU(Kl + l:NK> = FFIK) + UDI UNPQ) 
40 CONTINUE 

C EVALUATE FUNCTION AT ALL POINTS 


GLENS091 

GLENS092 

GLENS093 

GLENS094 
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APPENDIX B - Concluded 


c 

c 


NOTE CHANGE IN CALLING SEQUENCE FOR FUNCTION EVALUATION SUBROUTINE GLENS095 


CALL FXS(UUfFOFX.MAXI.NT) GLENS0V6 

CREATE VECTOR OF HEIGHTS OLENS097 

00 50 K*l,NO GLENS098 

KI=1K-11*NP0 GLENS099 

NK « K4' NPQ GLENSIOC 

UUU1«-1:NK) « RiJUNJ) GLENSIOI 

50 CONTINUE GLENSi02 

THESE STATEMENTS CAN BE REPLACED BY DOT PRODUCT FUNCTION IF AVAILABLE6LENS103 
SUCH AS GLENS104 

SUMin > DOTP (UUIlsNT) ,FOFX(lsNT. I) 6LENSI05 

GLENS106 


■** GLENS107 

00 80 loI.NFX 6LENS108 

SUMinsO.O GLENS109 

DO 70 K=liNT GLENSllO 

70 SUMin- SUMtII «'UUIK}* FOFXIK*!) GLENSlll 

SUMCIi > OELT * SUM(I) GLENSU2 

80 CONTINUE GLENSU3 

** GLENSlll 

RETURN GLENS115 

END 6LENSU6 
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